﻿using System;
using System.Collections.Generic;
using System.Text;

namespace GeoFly
{
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++
+														+
+				马斯京根－康吉法汇流程序				+
+														+
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+														+
+				Written by: ZhangDong					+
+				Date:	 2005-04-22						+
+														+
+ 计算步骤：1、设置汇流计算参数 SetRoutingPara();		+
+			2、计算汇流出流流量 RoutingOutQ();			+
+														+
+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
    /// <summary>
    /// 马斯京根－康吉法汇流程序
    /// </summary>
    public class MuskingumCunge
    {

        public MuskingumCunge()
        {
        }
        /// <summary>
        /// 设置汇流参数
        /// </summary>
        /// <param name="Vflow">流速</param>
        /// <param name="Slope">栅格坡度</param>
        /// <param name="deltaX">河长离散值</param>
        /// <param name="deltaT">时间离散值</param>
        /// <param name="riverType">河道类型</param>
        /// <param name="Q11">差分参数</param>
        /// <param name="Q12">差分参数</param>
        /// <param name="Q21">差分参数</param>
        /// <param name="dQ">当前时段的侧向单宽入流</param>
        /// <param name="dB"></param>
        /// <param name="x">从河源处到栅格处河道的长度</param>
        /// <param name="k">汇流时段</param>

        public void SetRoutingPara(double Vflow, double Slope, double deltaX, double deltaT, RiverType riverType,
                                  double Q11, double Q12, double Q21, double dQ, double dB, double x, double k)
        {
            m_dVflow = Vflow;
            m_dSlp = Slope;
            m_iRriverType = riverType;
            m_deltaT = deltaT;
            m_deltaX = deltaX;
            m_Q11 = Q11;
            m_Q12 = Q12;
            m_Q21 = Q21;
            m_dQ = dQ;
            m_dB = dB;
            m_X = x;
            m_K = k / 60;
        }
        /// <summary>
        /// 按Cunge法的有限差值法计算汇流流量输出
        /// </summary>
        /// <returns></returns>
        public double RoutingOutQ()
        {
            m_Q22 = 0;
            //下面计算系数
            double denominator = 1;
            denominator = m_K * (1 - m_X) + 0.5 * m_deltaT;
            C1 = (m_X * m_K + 0.5 * m_deltaT) / denominator;
            C2 = (0.5 * m_deltaT - m_X * m_K) / denominator;
            C3 = (m_K * (1 - m_X) - 0.5 * m_deltaT) / denominator;
            C4 = m_deltaT * m_deltaX / (denominator);
            //计算Cunge法的有限差分值
            m_Q22 = C1 * m_Q11 + C2 * m_Q12 + C3 * m_Q21 + C4 * m_dQ;

            //if (m_Q22 < 0)
            //    throw new Exception("汇流结果为负");
            return m_Q22;
        }

        /// <summary>
        /// 计算Cunge法的系数
        ///         KX + 0.5dt              0.5dt-KX                 K(1-X) - 0.5dt
        /// C0=  ----------------     C1=  ----------------     C2= ----------------    C4=...
        ///       K(1-X) + 0.5dt            K(1-X)+0.5dt             K(1-X) + 0.5dt
        /// </summary>
        public void CalcMuskRoutingCoeff()
        {
            C1 = C2 = C3 = 0.3333;
            m_dKWV = KinematicWaveV(m_dVflow);
            m_dDiffC = DiffusiveCoeff();

            /*	if(m_dKWV <= 0.)	
                {
                    m_X = g_ModelInputPara1.MCTimeWeight;
                    m_K = g_ModelInputPara1.MCGridRTravelTime / 60.;
                }
                else
                {
                    m_X = 0.5 - m_dDiffC / (m_dKWV * m_deltaX);	//时间差分权重
                    m_K = m_deltaX / m_dKWV;
                }
	
                if(m_K > m_deltaT)	m_K = m_deltaT;
            */
            double denominator = 1;
            denominator = m_K * (1 - m_X) + 0.5 * m_deltaT;
            C1 = (m_X * m_K + 0.5 * m_deltaT) / denominator;
            C2 = (0.5 * m_deltaT - m_X * m_K) / denominator;
            C3 = (m_K * (1 - m_X) - 0.5 * m_deltaT) / denominator;
            C4 = m_deltaT * m_deltaX / (denominator);
        }
        /// <summary>
        /// 计算扩散系数
        /// </summary>
        /// <returns></returns>
        double DiffusiveCoeff()
        {
            double dDC = 1;
            if (m_dSlp == 0)
                m_dSlp = 0.0001;
            dDC = (m_Q12 + m_dQ * m_dB) / (2 * m_dSlp * m_dB);
            return dDC;
        }
        double KinematicWaveV(double vflow)
        {
            double dKWVC = GetKWVCoeff(m_iRriverType);
            double dKWV = vflow * dKWVC;
            return dKWV;
        }

        double GetKWVCoeff(RiverType rivertype)
        {
            double dcoeff = 1;
            switch (rivertype)
            {
                case RiverType.M_RIVER_SECTION_TRIANGLE:
                    dcoeff = 1.33;
                    break;
                case RiverType.M_RIVER_SECTION_RECTANGLE:
                    dcoeff = 1.67;
                    break;
                case RiverType.M_RIVER_SECTION_PARABOLA:
                    dcoeff = 1.44;
                    break;
                case RiverType.M_RIVER_SECTION_HILLSIDE:
                    dcoeff = 3.0;
                    break;
                case RiverType.C_RIVER_SECTION_TRIANGLE:
                    dcoeff = 1.25;
                    break;
                case RiverType.C_RIVER_SECTION_RECTANGLE:
                    dcoeff = 1.50;
                    break;
                case RiverType.C_RIVER_SECTION_PARABOLA:
                    dcoeff = 1.33;
                    break;
                case RiverType.C_RIVER_SECTION_HILLSIDE:
                    dcoeff = 2.50;
                    break;
                default:
                    dcoeff = 1.5;
                    break;
            }
            return dcoeff;
        }



        double m_Q22;
        double m_dKWV;
        double m_dDiffC;
        double m_X;
        double m_K;


        double m_deltaT;
        double m_deltaX;

        double m_Q11;
        double m_Q12;
        double m_Q21;

        RiverType m_iRriverType;
        double m_dSlp;
        double m_dVflow;
        double m_dQ;
        double m_dB;
        /// <summary>
        /// 系数1
        /// </summary>
        double C1;
        /// <summary>
        /// 系数2
        /// </summary>
        double C2;
        /// <summary>
        /// 系数3
        /// </summary>
        double C3;
        /// <summary>
        /// 系数4
        /// </summary>
        double C4;
    }
}